{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1/6 + 5*sqrt(5)/6\n",
      "1/2\n"
     ]
    }
   ],
   "source": [
    "#3.1\n",
    "\n",
    "#求定积分符号解答，可以使用sympy中的integrate函数\n",
    "#(1)\n",
    "import scipy.optimize\n",
    "import sympy as sym\n",
    "x = sym.symbols('x')\n",
    "num = sym.integrate(sym.sqrt(1+4*x),(x,0,1))\n",
    "print(num)\n",
    "#(2)\n",
    "num = sym.integrate(sym.exp(-x)*sym.sin(x),(x,0,sym.oo))\n",
    "print(num)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "调用sympy求得的符号解为：[4/3 + (-1/2 - sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3) - 2/(9*(-1/2 - sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3)), 4/3 - 2/(9*(-1/2 + sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3), -2/(9*(64/27 + 2*sqrt(114)/9)**(1/3)) + 4/3 + (64/27 + 2*sqrt(114)/9)**(1/3)]\n",
      "{-2/(9*(64/27 + 2*sqrt(114)/9)**(1/3)) + 4/3 + (64/27 + 2*sqrt(114)/9)**(1/3): 1, 4/3 - 2/(9*(-1/2 + sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3): 1, 4/3 + (-1/2 - sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3) - 2/(9*(-1/2 - sqrt(3)*I/2)*(64/27 + 2*sqrt(114)/9)**(1/3)): 1}\n",
      "********************\n",
      "调用sympy求得的符号数值解为：[2.88123940107640, 0.559380299461801 - 1.56961032188996*I, 0.559380299461801 + 1.56961032188996*I]\n",
      "{2.88123940107640: 1, 0.559380299461801 - 1.56961032188996*I: 1, 0.559380299461801 + 1.56961032188996*I: 1}\n"
     ]
    }
   ],
   "source": [
    "#3.2\n",
    "\n",
    "# 注意如果要取出复数部分，只保留实数部分可以使用 solvest函数\n",
    "#另外求解器是会根据输入进行判别输出为什么类型的数据，比如输入的是整数，则输出完整的符号；输入的是浮点数，输出则为浮点数\n",
    "#scipy貌似对求解复数问题不适用\n",
    "\n",
    "#求普通方程的符号解和数值解\n",
    "# import sympy as sym\n",
    "# import numpy as np\n",
    "# x = sym.symbols('x',real = True)\n",
    "# expr1 = x**3-4*x**2+6*x-8#构造等式，然后再来解一元三次方程\n",
    "# num = sym.solveset(expr1,x)#求解根\n",
    "# print('调用sympy求得的符号解为：{}'.format(num))\n",
    "# print(sym.roots(expr1,x))#查看重根情况\n",
    "# import scipy as sci\n",
    "# f = lambda x : x**3-4*x**2+6*x-8\n",
    "# print(sci.optimize.fsolve(f,0))\n",
    "import sympy as sym\n",
    "import numpy as np\n",
    "x = sym.symbols('x')\n",
    "expr1 = x**3-4*x**2+6*x-8#构造等式，然后再来解一元三次方程\n",
    "num = sym.solve(expr1,x)#求解根\n",
    "print('调用sympy求得的符号解为：{}'.format(num))\n",
    "print(sym.roots(expr1,x))#查看重根情况\n",
    "print('*'*20)\n",
    "expr2 = x**3-4.*x**2+6.*x-8.#构造等式，然后再来解一元三次方程\n",
    "num2 = sym.solve(expr2,x)#求解根\n",
    "print('调用sympy求得的符号数值解为：{}'.format(num2))\n",
    "print(sym.roots(expr2,x))#查看重根情况"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "方程的符号解为： [(1/3 - sqrt(46)/3, sqrt(46)/9 + 17/9), (1/3 + sqrt(46)/3, 17/9 - sqrt(46)/9)]\n",
      "方程的符号化简数值解为： [(-1.92744332770842, 2.64248110923614), (2.59410999437509, 1.13529666854164)]\n",
      "Matrix([[-1.92744332770842], [2.64248110923614]])\n",
      "Matrix([[2.59410999437509], [1.13529666854164]])\n"
     ]
    }
   ],
   "source": [
    "#3.3\n",
    "\n",
    "#求方程组的符号解和数值解\n",
    "x,y = sym.symbols('x y')\n",
    "s = sym.solve([x**2-y-x-3,x+3*y-6],[x,y])\n",
    "print('方程的符号解为：',s)\n",
    "s = sym.solve([x**2-y-x-3.,x+3.*y-6.],[x,y])\n",
    "print('方程的符号化简数值解为：',s)\n",
    "s= sym.nsolve((x**2-y-x-3,x+3*y-6),(x,y),(-1,2))\n",
    "print(s)\n",
    "s= sym.nsolve((x**2-y-x-3,x+3*y-6),(x,y),(1,2))\n",
    "print(s)\n",
    "#发现一个很奇特的现象，如果你想求方程组的符号解，原方程输入即可\n",
    "#要求数值解，直接将数字变为浮点型恐怕更简单易用\n",
    "#另外，使用nsolve解函数时，提供的初值对于函数的输出会有影响，并且无法知晓函数到底有多少解\n",
    "#百度了很久，也没有合适的解答,这先存疑》》》》"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "边值问题的解为：Eq(y(x), -x*cos(2*x)/3 + (6*cos(4) - 4*sin(4) - 9*cos(2) + 27)*sin(x)/(9*sin(2)) + 4*sin(2*x)/9 + cos(x))\n"
     ]
    }
   ],
   "source": [
    "#3.4\n",
    "\n",
    "#求微分方程的符号解\n",
    "x = sym.symbols('x')\n",
    "y = sym.symbols('y',cls = sym.Function)#定义y为函数\n",
    "eq1 = sym.diff(y(x),x,2) + y(x) -x*sym.cos(2*x)\n",
    "print('边值问题的解为：{}'.format(sym.dsolve(eq1,y(x),ics={y(0):1,y(2):3})))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[1, 2], [3, 4], [5, 6]])\n",
      "Matrix([[1, 1], [2, 2], [3, 4]])\n",
      "Matrix([[1, 2, 1, 1], [3, 4, 2, 2], [5, 6, 3, 4], [2, 6, 3, 2]]) (4, 4)\n",
      "A的行列式为： 0\n"
     ]
    }
   ],
   "source": [
    "#3.5\n",
    "\n",
    "#求行列式的值\n",
    "A1 =sym.Matrix(np.arange(1,7).reshape(3,2))\n",
    "print(A1)\n",
    "A2 = sym.Matrix([[1,1], [2,2],[3,4]])\n",
    "print(A2)\n",
    "A3 = sym.Matrix([[2,6]])\n",
    "A4 = sym.Matrix([[3,2]])\n",
    "B1 = A1.row_join(A2)\n",
    "B2 = A3.row_join(A4)\n",
    "A = B1.col_join(B2)\n",
    "print(A,A.shape)\n",
    "print('A的行列式为：',sym.det(A))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A的基础解系为： [Matrix([\n",
      "[   0],\n",
      "[-1/2],\n",
      "[   1],\n",
      "[   0]])]\n",
      "(3, 4)\n",
      "(3, 1)\n",
      "增广矩阵的最简形为: (Matrix([\n",
      "[1, 1/2, -1/2, 0, 1/2],\n",
      "[0,   0,    0, 1,   0],\n",
      "[0,   0,    0, 0,   0]]), (0, 3))\n"
     ]
    }
   ],
   "source": [
    "#3.6\n",
    "\n",
    "#求解线性方程组,可以使用sumpy，也可以使用numpy\n",
    "#(1)其次线性方程+\n",
    "A1 = sym.Matrix([[1,2,1,-1],[3,6,-1,-3],[5,10,1,-5]])\n",
    "print('A的基础解系为：',A.nullspace())\n",
    "#（2）非齐次线性方程\n",
    "A2 = sym.Matrix([[2,1,-1,1],[4,2,-2,1],[2,1,-1,-1]])\n",
    "print(A2.shape)\n",
    "B = sym.Matrix([[1],[2],[1]])\n",
    "print(B.shape)\n",
    "C = A2.row_join(B)#构造增广矩阵\n",
    "print('增广矩阵的最简形为:',C.rref())#由最简矩阵可知，第1列和第4列为最简坐标向量\n",
    "#通过等价方程组，即可以解出最终的通解"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "系数矩阵的秩为： 2\n",
      "增广矩阵的秩为： 3\n",
      "系数矩阵的秩为： 2\n",
      "增广矩阵的秩为： 2\n",
      "最小范数解为： [[ 0.33333333]\n",
      " [ 1.33333333]\n",
      " [-0.66666667]]\n"
     ]
    }
   ],
   "source": [
    "#3.7\n",
    "\n",
    "#选择numpy\n",
    "#(1)非齐次线性方程组\n",
    "import numpy as np\n",
    "import numpy.linalg as LA\n",
    "A = np.array([[4,2,-1],[3,-1,2],[11,3,0]])\n",
    "B = np.array([[2,10,8]]).reshape(3,1)\n",
    "print('系数矩阵的秩为：',LA.matrix_rank(A))\n",
    "print('增广矩阵的秩为：',LA.matrix_rank(np.c_[A,B]))#说明系数矩阵的秩小于增广矩阵的秩，方程无解\n",
    "#（2）题目中有三个变量，但是有四个方程组，这种方程组称为矛盾方程组\n",
    "C = np.array([[2,3,1],[1,-2,4],[3,8,-2],[4,-1,9]])\n",
    "D = np.array([[4,-5,13,-6]]).reshape(4,1)\n",
    "print('系数矩阵的秩为：',LA.matrix_rank(C))\n",
    "print('增广矩阵的秩为：',LA.matrix_rank(np.c_[C,D]))#可见第二问有最小二乘解\n",
    "print('最小二乘解为：',LA.pinv(C).dot(D))#最小二乘解与最小范数解为同样的计算方式"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "特征值为： [ 2. 11.  2.]\n",
      "特征向量为： [[-0.74535599  0.66666667 -0.34869867]\n",
      " [ 0.2981424   0.33333333 -0.65103258]\n",
      " [ 0.59628479  0.66666667  0.67421496]]\n"
     ]
    }
   ],
   "source": [
    "#3.8\n",
    "\n",
    "#求特征值和特征向量\n",
    "A = np.array([[6,2,4],[2,3,2],[4,2,6]])\n",
    "values,vectors = LA.eig(A)\n",
    "print('特征值为：',values)\n",
    "print('特征向量为：',vectors )"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-a**2 + 2*a*b - 2*a - b**2 + 2*b - 1\n",
      "2*a*b + 2*b\n",
      "a**2 + 2*a*b + 2*a + b**2 + 2*b + 1\n",
      "a,b的值更新为 [(-1, 0)]\n",
      "A矩阵更新为: Matrix([[1, 0, 1], [0, 1, 0], [1, 0, 1]])\n",
      "特征值为： [2. 0. 1.]\n",
      "特征向量为： [[ 0.40824829 -0.40824829  0.        ]\n",
      " [ 0.          0.          0.57735027]\n",
      " [ 0.40824829  0.40824829  0.        ]]\n"
     ]
    }
   ],
   "source": [
    "#3.9\n",
    "\n",
    "#求参数a，b以及对应的正交变换矩阵\n",
    "'''\n",
    "二次型f(x，y，z)=ax**2+by**2+cz**2+dxy+exz+fyz，用矩阵表示的时候，矩阵的元素与二次型系数的对应关系为：A11=a，A22=b，A33=c，A12=A21=d/2，A13=A31=e/2，A23=A32=f/2\n",
    "由题可知，特征值为0,1,2\n",
    "利用两个性质：\n",
    "        1：设有N阶矩阵A，那么矩阵A的迹（用tr（A）表示）就等于A的特征值的总和\n",
    "        2：根据特征值为1和2，确定行列式A的值为0\n",
    "'''\n",
    "import numpy as num\n",
    "a,b = sym.symbols('a b')#real = True,constant = True\n",
    "A = sym.Matrix([[1,a+1,1],[a+1,1,b],[1,b,1]])\n",
    "print(A.det())#计算特征值为0时行列式的值\n",
    "print((A-sym.eye(3)).det())#矩阵的特征值为1行列式的值\n",
    "print((A-2*sym.eye(3)).det())#矩阵的特征值为2行列式的值\n",
    "exp1 = -a**2 + 2*a*b - 2*a - b**2 + 2*b - 1\n",
    "exp2 = 2*a*b + 2*b\n",
    "exp3 = a**2 + 2*a*b + 2*a + b**2 + 2*b + 1\n",
    "print('a,b的值更新为',sym.solve([exp1,exp2,exp3],[a,b])) #求解出a的值为-1，b的值为0\n",
    "print('A矩阵更新为:',sym.Matrix([[1,0,1],[0,1,0],[1,0,1]]))\n",
    "A = np.array([[1,0,1],[0,1,0],[1,0,1]])\n",
    "values,vectors = LA.eig(A)\n",
    "print('特征值为：',values)\n",
    "print('特征向量为：',vectors/num.linalg.norm(vectors) )"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cos(1)\n",
      "   2                \n",
      "  x ⋅sin(1)         \n",
      "- ───────── + cos(1)\n",
      "      2             \n",
      "                          2                \n",
      " 4 ⎛  cos(1)   sin(1)⎞   x ⋅sin(1)         \n",
      "x ⋅⎜- ────── + ──────⎟ - ───────── + cos(1)\n",
      "   ⎝    8        8   ⎠       2             \n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6w0lEQVR4nO3deViU5dcH8O8Igpr7ghq4IaAIoimilVkq4q5ZZqaZa2pqWpaWW2m5ZVaulaQmmWlZ5lroTzP3N8RcMiv3BXDBFUXZ7/ePE6ixyDIz9zMz3891cSkwzHOGZc7c2zkmpRSIiIiMppDuAIiIiLLCBEVERIbEBEVERIbEBEVERIbEBEVERIbk/IDPc4sfUR60adMG4eHhusMgsjWmrD7IERSRGV2+fFl3CER2gwmKiIgMiQmKiIgMiQmKiIgMiQmKiIgMiQmKiIgMiQmKiIgMiQmKiIgMiQmKiIgMiQmKiIgM6UGljogoK6mpQGwscO4cEBcHJCYCCQnA1avA118DRYsCpUoBVasC1asDLi66IyayOaYHdNQ1Ty2+114DDhwwy10RmcvVq1dx7PhxQClUrlwZVatWve/zl86exc0zZ1DaZEKRtDS4OjnBWSkgLQ3I5u8mEEBkdhc0mYBChYDChQFXV+Chh4ASJYAyZeR9IltSvz4wa5a57i3LWnwcQZFDUkrh6LFjqFevHlxdXXFg715UTEmBa1wccPs2kJwMNwBu935RaqokGBcXSTLpb05OknxMJuDSJaByZUliqalAUhKQknL/W0KCvN24cfe+TSa53xIlgAoVgPLl5VpEDsw6Iygig9mzbRv2DBmCkW5uwP79UDdu3H0J5+oKVKyIqJIlcaBoUXSYOxcIDJRE9ACBgYGIjMx2DHXX1avA9u3Arl3Avn3A0aPAhQuS1ABJWJUqAU2aAP36Ae3by8eI7FOWv9xMUOQ4rl4FZs8Gli+HOn4cpvTf/aJFcbliRUSUL492338PVKsGAFiyZAnGjBmDChUqwMfHB5988gmqVKmS6W5DQ0MRGhoKAIiNjcWZM2fyF59SwG+/AV9+CWzdCpw6JSMuAHB2limVYcOAF1/MVbIksiFMUOSAEhOBjz4CFi8GTpzI+HB8hQrYXaECWq1eDXh7Y+nSpfjtt98wb968jNtcuXIFxYsXh6urKxYsWIBvv/0Wv/zyS46Xy/UIKjeUkhHWZ58B//ufbMoAJFk1bAi8/jrw3HOcCiR7wH5Q5EB+/RV44gmgWDFg3DhJTl5ewHvvATdu4NCaNZjp4QF4ewMAoqKi4O7uft9dlCtXDq7/bl4YMGAA9u3bZ93HYDIBTZsCy5bJ2taJE0CfPkDJkjLS6t4dKF4c6NoVOH7curERWQETFNmPpCRg+nSgYkWgeXNg507ZIffGGzK9d+wYMGECULIkGjVqhGPHjuHUqVNISkrCihUr0KlTp/vu7vz58xn/X7t2LXx9fa39iO7n6SnTf1euAIcOAc88I6OsH36QROvjAyxdencdi8jGMUGR7Tt/HnjhBRlNjBkjU2FNmgDbtgGXLwMzZ0qiuoezszPmzZuH1q1bw9fXF926dYOfnx/eeecdrF27FgAwZ84c+Pn5oV69epgzZw6WLFmi4cFlo25dSUy3bwMrVgC1a0sCfukleawTJgB37uiOkqhAuAZFtuvYMeCVV4BffpGRRNGiQI8eMooqX15LSGZdg8qrf/6RdamNG2Wbe+HCQK9ewIcfAmXL6omJKHe4BkV24sgR4KmnZEpryxZZk/ngA6nosHChtuSkXa1awE8/yQjy5Zdl88TixYCbm4ysLl7UHSFRnjBBke345x9ZW/Lzk+m7UqVk2/jly8Do0bK7jWS0FBoKXL8OjBwpiWrpUsDDAxgxArh2TXeERLnCBEXGFxMjO9Z8fWV3XsmSdxPT8OFMTNkpUkS22J8/D/TuLZsn5swBHn4YmDgRiI/XHSFRjpigyLhu3JAdeNWqAd9+K4dTx4yRigtMTLlXrhywZIlsRW/dWsosTZokJZnmzLl7GJjIYJigyHjS0oBFi6QS+McfyxNo167A2bPA1KmyGYLyztMTCA8H/u//gHr1gJs3ZcrP11c2mhAZDBMUGcv//R/QoAEwYIBsevDykkOpK1fKK34quMaNgf37ge++k3W848eBli2BZ58F8lumicgCmKDIGC5ckHWSRx+VQ6iFC8to6cgRIChId3T2x2SSMkknTsgZMgD48UfZCfjuu3K+ikgzJijSKy1Nas2lV0EApETRkSOy3lS4sN747F25csA33wDr1sl29MREKQfl6yvTgUQaMUGRPn/9BTRrBgwZIq/YixaVc0y//ipTe2Q9HToAf/8N9O8v71+4ALRtC/TsebdILZGVMUGR9SUmyi6y+vWBvXvlY40aydRe//7se6RL6dLyAmHjRqlnCNwtoxQWlm0XYSJLYYIi69qzRzZBTJwo28RTUuT/O3YANWvqjo4AICQE+PNPWRNM7wzcp498/ORJ3dGRA2GCIutISADefBN4/HE5OFqokLxK37VLFuV5pslYSpSQs1OLF8vPrmRJYPduICBAqlRwNEVWwARFlrdvnzTY++gjGSVduyYtzH//XaqOk3H17QtEREj7+Tt35N9Bg+TnFxOjOzqyc0xQZDnJybIjrEkTKUtUtapMEU2ZAqxeLWseZHwBAUBkJPD887ItvU4d2cji7y8VPogshAmKLOPoUeCxx2T67tFHZZfe7duyAD92LNuU25oSJWQ7+qefysHeMmWAKlWkRmL37tJEkcjM+CxB5qWUrFs88oiMll54QTZA1KkjU3rBwbojpPwymaT/1q5d8v6pU7KR4ocfZJS1dave+MjuMEGR+Vy/Lq+m+/eXbeMtWgDLl0uS2rZNXnGT7QsMlPJTnp7A118Db70lI6yWLWV0nJysO0KyE0xQZB47d0oB0lWrgPHjZfv499/Leadly6T1A9kPDw8ZGYeEyJpi27ayoWLaNKBpU1mrIiogJigqmNRUOcf05JNSluibb+QtMlJGT++8w4O39qpECWDtWpn2mzVL2qN8/bU0lqxfX16YEBUAD59Q/p0/L6Vwtm6VluJ9+0oBUpNJdnlxC7n9c3YG5s8HvL2ld1dUFLB5M/Daa8CLL8rvwZw5bJFC+cIERfmzZYskp7g4OdD58MNSz61CBWDTJnnCIsdgMgGvvw7UqAH06CEvUtaskTJJ06bJetV330nJJKI84BQf5U1qqqwrtWoFlC0rtfSKFJGDmzVrSrUBJifH9PTTwPbtUnmieXPpL/XzzzLSDgzklB/lGRMU5d7Fi9IyfOJEoFcvSU5bt8ouvSZNZKcemwo6tsBA2YZesqTs6itVCjhwQOovvvgi8PLLUpGCKBeYoCh3tm2The/du+WcU1gY8MknwKuvytTexo02VxkiPDwctWrVgpeXF6ZPn57p84mJiXj++efh5eWFxo0b4/Tp09YP0hZ5espIqkIFGWkfPy4t5ceMkWrpjz3GXX6UO0qpnN7I0aWlKTVjhlJOTkrVqqXUH3/IxyZMUApQqlcvpZKTdUeZZykpKcrT01OdOHFCJSYmqoCAAPXnn3/ed5v58+erQYMGKaWUWr58uerWrdsD77dhw4YWidcmxcQoVaeOUkWLKrVxo3xs3TqlSpdWqlQppVav1hoeGUqWOYgjKMrejRuyjjB6NNCli0zp+fnJYcz33wf69QO+/NImK5FHRETAy8sLnp6ecHFxQffu3bFmzZr7brNmzRr07t0bANC1a1ds2bIFilW8c69yZdnF5+MDdOwoXXs7dJCKIjVryprVW2/JmTmiLJhy+oNr06aNunz5slkuFBsbiwoVKpjlvmyJzT7uO3dkGiYxUQ5lpjewi4qStajy5YFq1bL9cqM/7mvXriEuLg7V/n0MV65cQXx8PKpWrZpxmz///BPe3t5wcXEBAPzxxx/w9fWF838ScmxsLNL/ThITE1G/fn3rPAiDyfZnnpICHDsmtRg9PaWOn1LA2bNSRLh4cfl44cLWD9oMjP67binmfNz79u3bqJRqk+kT2Q2tlJmn+Bx16sMmH/fSpTItU7myUjt2yMfS0pQaPlym9YYNk/dzYPTHvXLlStW/f/+M97/66is1dOjQ+27j5+enzp07l/G+p6enio2NzfF+ixUrZt5AbUiOP/Pr15V6/HGlChVS6uuv7348LEypIkWUevhhpXbvtnyQFmD033VLMfPj5hQfPUBSEjBsmOzQCwqSqZimTeXV7ujRcuDy9dflXxuvDuHu7o5z585lvB8VFQV3d/dsb5OSkoIbN26gXLlyVo3TbpQqBYSHS8WRXr2k2gggB7z37JGjCk8+KdXSOY1K/2KCIhEdDTz1lFQFePNNqQZQqZJ8buJEYOZMYOhQaTpo48kJABo1aoRjx47h1KlTSEpKwooVK9CpU6f7btOpUyeEhYUBAL7//nu0aNECJjt47NoULw5s2CCJ6KWXZE0KkN2hkZFS12/oUKmQfvu21lDJILIbWikzT/EtWLDAnHdnM2zicW/bplTFiko99JBS3313/+emTZNpvf79lUpNzfVd2sLj3rBhg/L29laenp5q8uTJSimlJkyYoNasWaOUUurOnTuqa9euqmbNmqpRo0bqxIkTD7zPqlWrWjRmI8v1zzwuTqlGjZRydVVqy5a7H09NVeq995QymZQKCFAqF99vI7CF33VLMPPjzjIH5bhJAgDH2vZMKWD2bBkx1awJ/Pij9G1KN2cOMGKElK/56ivAyUlfrDYiMDAQkZGRusMwvitXZMR+6pSUzWrc+O7nwsPld04pqT7Rrp22MMlqspya4BSfo4qPl1p6r78uW4AjIu5PTl99JcnpmWfkUC6TE5lTuXJSs7FSJWnV8ccfdz/Xpo1M+VWvLtvS33sPSEvTFirpwwTliI4dk9JE334LTJ0qHVFLlbr7+Q0b5IxTy5aymG2D55zIBlSuLGudxYrdrTiRztNTSia9+CLw7rtA587SEJMcilUT1IQJExAQEID69esjJCQEMTEx1ry8NqNGjULt2rUREBCALl264LrOP7R166Tb7fnzUshzzBig0D2/Brt3SzXq+vVlys/VNd+XWrlyJfz8/FCoUCGHmPYKDw/H4cOHsy2dZK/69esHNzc3+Pv75/2Lq1eXJJWaCgQHA/fsrESxYjJ6nzdPpv0CA+8faWl07tw5NG/eHHXq1IGfnx9mz56tOySrSUhIQFBQEOrVqwc/Pz+8++67lrtYdotTygKljm7cuJHx/9mzZ2eUkbF3GzduVMn/lgMaPXq0Gj16tPWDSElRavx42fDQsKFSp09nvs3hw0qVKaOUt7dSFy8W+JJHjhxRf//9t3ryySfV3r17C3x/RpZeOsnf3z/b0kn2atu2bWrfvn3Kz88v/3eyb59SJUsq5eur1LVrmT+/c6ecyytWTKlvvsn/dcwkJiZG7du3TymlVFxcnPL29naYn3daWpq6efOmUkqppKQkFRQUpPbs2VPQu9V/DqpkyZIZ/4+Pj3eYLbshISEZ1QeaNGmCqKgo6wZw9aq0w5g8Wabudu7MXAXi7FmpVF6kiKwNuLkV+LK+vr6oVatWge/HFqSXTnJ1dc22dJK9atasGcqWLVuwO2nQQLrzHj8u5bWSku7//OOPy7m8Bg1kA8WIEZlvY0WVK1dGgwYNAAAlSpSAr68voqOjtcVjTSaTCcWLFwcAJCcnIzk52WLP5dZJUBcvAteuAQDGjRuHKlWqYNmyZXjvvfescnkjWbx4Mdq2bWu9C+7fDzRsKG0xFiyQatJFitx/m7g4SWA3b8pUSvXq1ovPTkRHR6NKlSoZ73t4eDjME5bZPPkksGiRVD5/5ZXMB3YrVZLPjRghO0ybNwcMsExw+vRp7N+/H43v3Ylo51JPnEBQQADc3NzQqlUriz12syeo4OBg+Pv7Z7wF1qmDSx4eiOrcGQAwZcoUnDt3Dj179sS8efPMfXlt/vu409/ufRU9ZcoUODs7o2fPntYJatEi4NFHpRbajh3AwIGZD9kmJ8ua099/y2aJgIA8XSI3j5so13r1At55R1q6ZLWOV7gwMGsWsHw5cPCgjKi2bbN6mOlu3bqFZ599FrNmzbpvhsiuJSbC6emnEVGxIqKiohAREYHDhw9b5lrZzf0pc65Bvf++rH18+23Gh86cOVOwOWsb8+WXX6omTZqo+Ph4y1/s9m2l+vWT73lwsFKXLmV9u7Q0pQYNktstXGixcBxhDWr37t0qJCQkoz7Z1KlT1dSpUzVHZT2nTp0y399zWppSPXpkes7I5PBhaQHj5KTUhx8+sD6kuSUlJamQkBD10UcfWfW62o0bJz+bDRuUUkpNmjRJffjhhwW91yxzkHUSVHKyUkFBKqV0aekRo5SaM2eOevbZZ812CSP7+eefla+vr7qUXaIwpxMnlHrkEfnRjh8vmyOy8+GHcrsxYywakiMkqOTkZFWjRo37NkkcPnxYd1hWY9YEpZRSCQlKNW0q1SZyKiJ744ZSzz4rv8fPPCPvW0FaWprq1auXGjFihFWuZxh796o0JyeV0KOHUkqp27dvq6ZNm6p169YV9J41JiillPr7b5Xg5KS2FS+u6vr7qw4dOqioqCizXsKoatasqTw8PFS9evVUvXr1LLd7Mb0ZXOnSSq1fn/Ntv/9efvzduuWphFFerFq1Srm7uysXFxfl5uamQkJCLHIdo9iwYYNydXW9r3SSI+jevbuqVKmScnZ2Vu7u7mqhuUbjsbFK1aypVPnyOZc9SktTauZMGUl5eyu1f795rp+DHTt2KACqbt26GX/XG/4dUdithASl/PxUUoUKqqm/v6pbt67y8/NTkyZNMse9a05QSik1Z45cMjTU7Hft0JKSlHrzTfnePvLIg2uY/fabtDh49FGZDiSzcdTWCxbzzz9y9KF2baWuXs35ttu3S9sOV1elPv/c6lN+dm/s2Pum9szMAAkqNVWpFi2kKKmNFII0vDNnJNEASg0ZotSdOznfPjpazpPUqJH92hTlGxOUBWzbplThwko1b65UYmLOt714UamQEPl76NFDqX/P61ABbd8uI9Q+fSx1Bf3noFCokLQId3IC+vSR0+OUf+vXA488Ahw+DKxYIa0y/ruF/F4JCVJbLy4OWLMGcMAuoGSDmjWTHalbtwKDB+fcL8rNTSqkTJ4sfxMGqj5hsy5dArp3B2rUkOLSVmT9WnxVq8oZhh07gA8+sPrl7UJysjQQ7NhRvp/79gHPP5/z1ygFDBkC/PabFIKtW9c6sRKZQ/r28y+/BObOzfm2hQoB48ZJlfQbN6T55qJFbISYH6mp8r2/cgVYuRKw9lb67IZWyhJTfOnS0pR6/nlp/7xxo8UuY5eOH1eqcWOZwnjllQdP6aVLX/975x3LxufgOMVnQampSnXsqJSzs1I7duTuay5cUKplS/nd79496zJKlL333rPWvgEDrEHd69Ytpfz9lSpbVqmTJy16KbuQlqZUWJhSxYsrVapUzudD/mvLFpk/7tzZYjv2SDBBWdi1a0p5eSlVqVLGkZUHSklRavJkSWxVqij1668WDdFubNkig4gXX7TGhhMDrEHd66GHpFp2aqqsi7DFc/Zu3JDeTb17y5rToUNAt265+9rTp6VSRK1awNKl91cuJ7I1pUvL80ZcnPxe56Yen5OTTPnt3i1rtM2bA2+/rbWWn+GdPCnLBj4+wGefZa5AYyV6n628vKRj5sGD0veFmyYy27ULqFcP+O47WfjdulXWnXIjIUEKb6amyqaIEiUsGyuRNfj7y5rSrl3AqFG5/7pGjaQ25csvy/r3o49KiS+639Wr0sU4LU0K+P5bGFaL7IZWytJTfPeaNUtmGx3tVHZO7txRavRoGWJ7eiqVn3L2AwfK93XNGvPHR/eJiIhQdevWVQ0aNFC3bt1SderUUX/88YfusOzba6/J7/eyZXn/2tWrlSpXTqmiRZX69FOemUqXmKjUk08q5eIiW8utx2BrUP+V/ss2Y4ZVL2tIERHSFweQJBMXl/f7CAuTr3/7bfPHR1kaN26cqlixohoyZIhD1eHTJilJqWbNJMkcPJj3r4+JUap1a/k7addOqbNnzR+jLUlJkY0k+U36BWPwBJWSImV3AKVmz7bqpQ0jIUEKMTo5KeXurlR4eP7u5+BB+aNt3lzqIJJVJCYmqqJFi6qgoCCVklMNRDKf8+elekTNmvnboZeaKjtcixWTDUhz5uRcv9JepaYq9dJLOgcJBk9QSskroi5dJKx586x+ea22b1fKz08ee58++d8Oe/267HKqXFm22JLVxMTEKBcXF+Xr66tu3bqlOxzHsWuXVJro0CH/u1RPnbo7mmrcWClHmp5NTVVqwAB57O+/rysKG0hQSskcaOfOEtr8+VpCsKrY2LutMapWfXCR15ykpUlFZyen3J8TIbPp2LGjqlGjhpo8ebIaOnSo7nAcy7x58jc0bVr+7yMtTamvv5bitM7OMpuR23OGtiolRamXX1YZ3Q/0sZEEpZQkqU6dJLwJE+xzATM1ValFi2Sh1tlZqbfekrNhBTF/vnzPZs40T4x26sqVKyo4OFh5eXmp4OBgdTWbIqSFChXKqFTdsWPHHO8zLCxMPfPMM6phw4YqJSVFBQUFqS1btlgifMpKWposETg5yYiqIGJj7053+fjY77mp27eVevppeZxjx+p+nrWhBKWUrJ307393yispSWs4ZnXggPS6AeRfc0wnHDwoVZzbtuVh3AcYNWqUmvbvK+1p06ap0aNHZ3m7hx56KM/3zYO6Gl2/LkWQq1Z9cOXz3Ni4Ue4vvQqFPRW4vnxZqccfV8pkknU3/WwsQSklGX3iRAnziSdyf3LcqM6ckVdmJpOMnBYvNk8yuXVL2hFUqiTVnClHPj4+Kubf36WYmBjl4+OT5e2YoGzQb7/JjESXLuYZEdy6JVN9RYvKOteIETLCsmWRkUpVqyZbyb/7Tnc06WwwQaVbtkx22VSqJKX3bc2VK9KvydVV3kaPNs8rvHT9+0vS27zZfPdpx0qVKpXx/7S0tPvev5eTk5Nq2LChaty4sfrxxx+zvb8FCxaohg0bqoYNG6qqVauaN1jKu5kzzb+GHRUlGwkKFVKqZEmlpkxRKj7efPdvLV98Ic9BVapIMjcOG05QSsk0mI+P/IKMGmUbi5dxcUp98IF0uDWZZKrS3GctVqxQ1mjbbmtatmyp/Pz8Mr2tXr06U0IqXbp0lveR3vH5xIkTqlq1aur48eMPvC5HUAaQmirnmlxdzd9d988/766PP/ywUgsX2sZRjthYKdANKNWqlRFHgTaeoJSSJ/z06gi1ayu1e7fuiLIWFSWjpFKlVMYhwEOHzH+dkyfl1VyTJva1RmdhuZ3iu1fv3r3VypUrH3g7JiiDiI2VBFKrVsE3H2Vl+/a7XQW8vZVasMCYL5rT0pRavlypihVlinLyZKOe88oyB9lW5dASJYAFC4BNm4D4eOCxx4CXXgKionRHJv74Qxox1qgBzJwJtG4t/Zc2bDB//6XkZOCFF6SI4/LlQOHC5r1/O9apUyeEhYUBAMLCwtC5c+dMt7l27RoSExMBAJcvX8auXbtQp04dq8ZJBVC+PPD118DRo3mr15dbTzwB7NkD/PCD9EgaNAioXl3qZV64YP7r5cfvvwMtW8rzhLs7sHevFM11ctIdWe5ll7mUEUdQ94qLk2ktFxd5e/VV2YRgbbdvK/XNN3d7zhQrJrFYuoXI+PFyPeMsctqMy5cvqxYtWigvLy/VsmVLdeXKFaWUUnv37lX9+/dXSim1a9cu5e/vrwICApS/v79auHBhru6bIyiDefNN+TspyPnCB0lLU+qXX5Rq00au5ewsW943b9YzWtm/X6muXSWWcuVkLc6Yo6Z72cEUX1ZOnZLFSycneXv+eeljYsmt1omJSm3YoFTfvrK+BMiumClTZEOEpe3ZI2txvXtb/lqUJ0xQBpOQoFRAgFJubtbZ4fr330qNHKlUmTLyvFCpklLDhim1c6dln5MSEpRauVLKmwEy9T9unGy9tw1Z5iCTUjm2QbadHslnz0or6IULgevXgSpVpM9Uhw5A06bSBya/UlOBv/4Ctm+X6cVffgFu3pShfefOMq331FPW6bUUHy89oRITpS9UqVKWvyblWmBgICIjI3WHQfc6fBgIDARCQqTtjDV6G925A6xfD3z7rUzxJyQAbm5AixZAcLBMvVWrVrBYrl2T56N16+Qa169LK54hQ2TKsXRpcz0aa8jyG2E/CSrdnTvA6tXAihVAeLg0JStcGGjQQJ7Y/fzkh1i5MlChAuDiIp9XSn7g164Bly4Bx4/L/PU//wD79klCAmSeOSQE6NRJftFcXa37+IYOBT79VPpCPfWUda9ND8QEZVCzZgGvvy5r2AMHWvfaN29KX6XwcGDz5rtrVG5uQP36QEAA4Okpzy2VKskL3xIl5AVvYqIkt4sXgeho4MwZ4M8/gchI+VcpWW9r316amrZoYVtrTHc5SIK6V3w88OuvwI4d0k3z0CHpTptbZcpIJ9oGDYAmTWRThqentu6S2LgRaNNG/tA+/lhPDJQjJiiDSkuTTUu7d0vTQh8fPXEoJbMxW7fKC98DByTR5KW7b6VKktgee0xGYo0b22pSupcDJqj/UkpevURFAefPA5cvy264lBT5fOnSkpTKlZNuv+XKaQ33PlevSifRMmXkF7sgU5ZkMUxQBhYdLbtpvbwkUTk7645IpKXJ89Hp0zJ7c/OmvKWlyQyNiwtQsSLg4SFvZcrojtgSskxQBvkJWYnJJFN7lSvrjiTvhgwBYmNlrpnJiSjv3N2Bzz8Hnn8emDEDGDtWd0SiUCGJzd1ddySGY1vnoBzVihWy2DpxoqyjEVH+dOsGPPec/C0dPqw7GnoAx5ris0XR0TK1V7u2rKUZZVqCssQpPhtw6ZJslqpeXQ7b8m/KCLKc4uMIysiUku2iiYnAV1/xD4nIHNzcZCdsZKRM9ZFhMUEZ2fLlsuY0ZQrg7a07GiL78dxznOqzAZziM6pLl4A6dSQx7dxpD9tIHQKn+GwIp/qMhFN8NmX4cNlqumgRkxORJXCqz/CYoIxozRrZtTdhgoyiiMgyONVnaJziM5rr1yUpVaggr+zYRsOmcIrPBqVP9dWoIVN9nLHQgVN8NuHNN+UPZvFiJicia3BzA+bMkX5J8+bpjobuwQRlJFu2yJrTm28CDRvqjobIcXTvDrRtKw39zp7VHQ39iwnKKOLjgZdfll17776rOxoix2IyyYYJpaSsWM5LH2QlTFBGMX48cOqUjKCKFtUdDZHjSW/ZvmEDsHKl7mgITFDGsHcvMHu2vHJ74gnd0RA5rldflen14cOlNxxpxQSlW2oqMHiw9HiZNk13NESOzdkZ+OILacXz1lu6o3F4TFC6ffop8Pvv0vGzZEnd0RDRI48AI0dKotq2TXc0Do3noHQ6f16qlDdpIu2gdXXqJbPhOSg7cfu2dBEoXFg6cbu66o7I3vEclOGMHCmVyufPZ3IiMpJixaS54dGjwMyZuqNxWExQuvzvf9KIcOxYaUFNRMYSEgJ07SrdBE6f1h2NQ2KC0iEhQXbseXtzIVaDlStXws/PD4UKFcpxOi48PBy1atWCl5cXpk+fbsUIyTA+/lhasr/2mu5IHBITlA4ffAAcPy4bJDi3bXX+/v5YtWoVmjVrlu1tUlNTMXToUPz88884cuQIli9fjiNHjlgxSjKEKlWAd96RAs4bNuiOxuEwQVnbsWPA1KnACy8AwcG6o3FIvr6+qFWrVo63iYiIgJeXFzw9PeHi4oLu3btjzZo1VoqQDOW11wBfXzkbdeeO7mgcChOUNSkFDB0KFCkiUwdkWNHR0ahSpUrG+x4eHoiOjs7ytqGhoQgMDERgYCBiY2OtFSJZi4uLFJE9eZJ9o6yMCcqaVq6UzRFTpsjBXLKY4OBg+Pv7Z3qzxCho4MCBiIyMRGRkJCpUqGD2+ycDaNFCCspOmwacOKE7GofBHsfWEh8PvPGGHAJ85RXd0di9zZs3F+jr3d3dce7cuYz3o6Ki4O7uXtCwyJbNnAmsXy9TfevX82iIFXAEZS0ffABERUnfGTZEM7xGjRrh2LFjOHXqFJKSkrBixQp06tRJd1ikk7s7MGkS8NNPwLp1uqNxCExQ1nD6NPDhh7IxomlT3dE4vB9//BEeHh7Ys2cP2rdvj9atWwMAYmJi0K5dOwCAs7Mz5s2bh9atW8PX1xfdunWDn5+fzrDJCF59VbrvDh8u1SbIoljqyBq6dgV+/hn45x/Aw0N3NGRBLHXkALZtA556CpgwAXjvPd3R2AuWOtLil1+AH34AxoxhciKyB08+KRsmPvwQOHNGdzR2jQnKklJSgBEjgBo1pI07EdmHDz6QTRKsBGNRTFCW9PnnwOHDwEcfydknIrIPVasCo0YB334L7NypOxq7xTUoS7lyRWrtNWggZ5+4JdUhcA3KgcTHA7VqyZnGiAip2Uf5xTUoq5owAYiLk1buTE5E9uehh2Sqb98+ICxMdzR2iQnKEg4eBBYskIrl3JpMZL969JCGo2PHAjdv6o7G7jBBmZtSUlyyTBk51EdE9stkklmSCxekCDSZFROUua1dC/z6q5yPKFNGdzREZGlBQUCvXlIA+uRJ3dHYFW6SMKfkZMDfXxZLDx0CChfWHRFZGTdJOKjoaMDHB2jTRs49Ul5xk4TFLVgAHD0qJfmZnIgch7u7HMZftQrYsUN3NHaDIyhzuXEDqFkTCAgAtmzhzj0HxRGUA7t9W0ZRHh7Anj18DsibLL9ZbLdhLlOnAlevSkl+/mKa1aVLl7Br1y7ExMSgaNGi8Pf3R2BgIArx3AkZSbFiwPvvA/36Ad9/Dzz3nO6IbB5HUOZw+jRQuzbw/PM8D2FGW7duxfTp03H16lU88sgjcHNzQ0JCAo4ePYoTJ06ga9eueOONN1CyZEndoWbgCMrBpaYC9etLa/gjR6QbL+UGR1AWM3asjJomT9YdiV356aef8MUXX6Bq1aqZPpeSkoL169fjf//7H5599lkN0RFlwclJ1qDbtZNSZ8OH647IpnEEVVAREUDjxsC4cUxQxBEUyVnIVq2AAwekPXypUrojsgXcxWd2Skkbdzc3VjW2oF69euHGjRsZ758+fRotW7bUGBFRDkwmGUVduQJMn647GpvGBFUQ69ZJJeNJk4ASJXRHY7eaNm2Kxo0bZ0z5hYSE4LXXXtMdFlH2GjQAevYEZs0Czp3THY3N4hRffqWmAvXqAUlJwJ9/8tyThe3cuRPNmzdH+fLlsX//flSqVEl3SFniFB9lOHNGtp336AF8+aXuaIyOU3xmtWyZJKbJk5mcLGzp0qXo168fvvrqK/Tp0wft2rXDwYMHdYdFlLNq1WSTRFiYFJCmPOMIKj8SE2VbedmywN697ANjYU8//TRCQ0Ph5uYGAIiIiMCgQYOwf/9+zZFlxhEU3efaNTnAHxQEhIfrjsbIOIIym9BQOfs0dSqTkxWsXr06IzkBQFBQEH777TeNERHlUpkywPjxwMaNwObNuqOxOXx2zaubN+W0+FNPASEhuqOxa5MnT8bVq1ez/JyLiwt++eUXrF+/3spREeXR0KFAlSpyFCXnGSv6Dx7UzatZs4DYWGDaNJY0srC6deuiY8eOKFKkCBo0aIAKFSogISEBx44dw4EDBxAcHIyxY8fqDpMoZ66uwLvvAgMGSDuezp11R2QzuAaVF5cvy3xyixbAjz/qjsbu9erVC0uXLsWMGTPg5uaG8+fPo2jRovD19UWzZs1QtGjRfN3vypUrMXHiRPz111+IiIhAYGBglrerXr06SpQoAScnJzg7O+dqbYlrUJSllBSgTh1JVgcPcmkgM5Y6KrDp04Fbt1gxwkr27duHmJgYLFu2DFu3br3vc3fu3Ml3gvL398eqVaswaNCgB95269atKF++fL6uQ5TB2VmamL7wArBihWw9pwdigsqtc+eAefOkc6afn+5oHMLgwYPRsmVLnDx58r5RjlIKJpMJJ/PZvdTX19dcIRLlXrdusjTwzjtS6ZzHUx6I48zcmjRJFjgnTtQdicMYPnw4/vrrL/Tr1w8nT57MeDt16lS+k1NemEwmhISEoGHDhggNDc32dqGhoQgMDERgYCBiY2MtHhfZqEKFgClTpD7fkiW6o7EJXIPKjb//llHTq6/KJgkyvODgYFy4cCHTx6dMmYLO/y5SP/XUU5g5c2a2a1DR0dFwd3fHpUuX0KpVK8ydOxfNmjXL8bpcg6IcKQU89pjMyBw/DhQpojsio+AaVL5NmCDNyLhjzGZsNsOZE3d3dwCAm5sbunTpgoiIiAcmKKIcmUxyfrJFC+Czz4DXX9cdkaFxiu9BIiOlO2Z61XJyCPHx8bh582bG/zdt2gR/f3/NUZFdaN4caNlS1qP+/R2jrDFBPciECUC5csDIkbojITP58ccf4eHhgT179qB9+/Zo3bo1ACAmJgbt2rUDAFy8eBFNmzZFvXr1EBQUhPbt26NNmzY6wyZ7MmWKnKecPVt3JIbGNaic7Nkj88XTp7PfE+UK16Ao155+Gti6FTh1Sup6OjbW4suziROB8uWlVAkRkTm9/75M8X34oe5IDIsJKju7dgGbNsnIqXhx3dEQkb2pW1cO7s6eDWSx45SYoLL37ruyKeKVV3RHQkT2atIkaXo6Y4buSAyJCSor27cDW7bI6Omhh3RHQ0T2yssLePFF2XLOUVQmTFBZefddoFIlYPBg3ZEQkb0bPx5ITuZaVBaYoP7r11/l7e235XAuEZEleXkBPXtyFJUFJqh7KSWjp4cfBgYO1B0NETmK8eOBxESOov6DCepev/wi609jxgD5bOVARJRn3t5316IuXtQdjWEwQaVLHz25u0vnSyIia+IoKhMmqHSbN8vZp3HjWGGYiKzP21vWoj79lKOofzFBATJ6mjQJqFIF6NdPdzRE5KjSR1EzZ+qOxBCYoABg2zYZPb31FuDqqjsaInJUPj4yipo/H7h0SXc02jFBAVITq1IloH9/3ZEQkaPjWlQGJqjdu2X33qhRXHsiIv18fIAePWQtysFHUUxQkydLxfJBg3RHQkQkxo8H7twBZs3SHYlWjp2gIiOBn3+WZoSsuUdERlGrFtC1q6xFXb+uOxptHDtBTZkClC7Nfk9EZDxjxwJxccC8eboj0cZxE9ShQ8Dq1cCIEUDJkrqjISK6X/36QPv2Ms0XH687Gi0cN0FNnQqUKAEMH647EiKirI0bB1y5AoSG6o5EC8dMUP/8A3z3nUztlS2rOxoioqw9+ijQvLkc3E1M1B2N1Tlmgpo+XbaUjxypOxIiopyNHQvExABLluiOxOocL0GdPQt8/TXw8stAhQq6oyEiylnLlkBQEPDBB0BKiu5orMrxEtRHH8m/b7yhNw4iotwwmWQt6tQpYPly3dFYlWMlqNhY4IsvpO9K1aq6oyFNRo0ahdq1ayMgIABdunTB9WzOmYSHh6NWrVrw8vLC9OnTrRsk0b06dADq1gWmTQPS0nRHYzWOlaDmzAESEoDRo3VHQhq1atUKhw8fxqFDh+Dj44Np06Zluk1qaiqGDh2Kn3/+GUeOHMHy5ctx5MgRDdESAShUSNai/vpLjsc4CMdJUDdvyoG3p58GfH11R0MahYSEwNnZGQDQpEkTREVFZbpNREQEvLy84OnpCRcXF3Tv3h1r1qyxdqhEdz33HODlJZu8lNIdjVU4ToJasEBKhowZozsSMpDFixejbdu2mT4eHR2NKlWqZLzv4eGB6OjoLO8jNDQUgYGBCAwMRGxsrMViJQfn5AS8+Sawdy/w66+6o7EKx0hQiYnAxx/LbphGjXRHQ1YQHBwMf3//TG/3joKmTJkCZ2dn9OzZs0DXGjhwICIjIxEZGYkK3BlKltS7N+DmBsyYoTsSq3DWHYBVhIUB588DS5fqjoSsZPPmzTl+fsmSJVi/fj22bNkCk8mU6fPu7u44d+5cxvtRUVFwd3c3e5xEeVKkiJRnGzcOOHgQqFdPd0QWZf8jqNRUafwVGAi0aKE7GjKA8PBwzJgxA2vXrkWxYsWyvE2jRo1w7NgxnDp1CklJSVixYgU6depk5UiJsvDKK0Dx4g4xirL/BLV6NXD8uLRzz+KVMjmeYcOG4ebNm2jVqhXq16+PwYMHAwBiYmLQrl07AICzszPmzZuH1q1bw9fXF926dYOfn5/OsIlEmTLAwIHAt98Cp0/rjsaiTCrn3SC2vVVEKalldfmy1N9zctIdEdm5wMBAREZG6g6D7F1UFFCjhoym5szRHY05ZDl6sO8R1M6dwG+/Sc09JicishceHkDPnsDChfIC3E7Zd4L68ENp596nj+5IiIjMa9QoaQs/f77uSCzGfhPUX38B69ZJS41sFsKJiGyWn5+UQJo7124bGtpvgvr4Y9mSyXbuRGSv3npLGhp++aXuSCzCPhPUhQvAV18BffuypQYR2a+mTYHHHpMuDXbYisM+E9TcuUByMvD667ojISKyrFGjZLv5qlW6IzE7+0tQt24Bn30GdOkCeHvrjoaIyLI6dpQish99ZHdFZO0vQS1aBFy7Jq8qiIjsnZMT8NprQEQEsHu37mjMyr4SVEoK8MknMi/bpInuaIiIrKNPH6kw8fHHuiMxK/tKUCtXAmfOcPRERI7loYeAwYOBH38ETpzQHY3Z2E+CUkpePdSqJWcDiIgcybBhgLMzMGuW7kjMxn4S1K5dQGSkzMUWsp+HRUSUKw8/DLzwArB4sazD2wH7eSb/5BOgbFngpZd0R0JEpMfIkcDt29JB3A7YR4I6dUraagwaxLJGROS46tWTzuFz5wJJSbqjKTD7SFBz5si0HssaEZGje+MNICZG+kXZONtPUDduyNmn558H2JKbiBxd69aAr69sGrPxg7u2n6AWLQJu3mRZIyIiQGaTRo4EDhwAtm7VHU2B2HaCSkmR6b0nngAaNtQdDRGRMbz4ohTKnj1bdyQFYtsJavVqOZjL0RMR0V1FisimsXXrbPrgrm0nqE8+AWrUADp10h0JEZGxvPKK1OmbN093JPlmuwkqvTDi8OHyQyAiorsefhh47jk5uHvzpu5o8sV2E9Ts2UCJEkC/frojISIyphEjgLg4ICxMdyT5YpsJ6vx54LvvJDmVLKk7GiIiY2rcGAgKks1kaWm6o8kz20xQCxYAqalSHJEoj0aNGoXatWsjICAAXbp0wfXr17O8XfXq1VG3bl3Ur18fgYGB1g2SyFxGjACOHQM2btQdSZ7ZXoJKSgI+/xxo21a6SBLlUatWrXD48GEcOnQIPj4+mDZtWra33bp1Kw4cOIDIyEgrRkhkRl27ApUr2+SWc9tLUCtXAhcvyuYIonwICQmBs7MzAKBJkyaIiorSHBGRBbm4yI6+jRuBv//WHU2e2F6CmjMH8PEBWrXSHQnZgcWLF6Nt27ZZfs5kMiEkJAQNGzZEaGhotvcRGhqKwMBABAYGIjY21lKhEuXfoEGSqObO1R1JnphUzrWajFXIKSJCFv3mzuX6E+UoODgYFy5cyPTxKVOmoHPnzhn/j4yMxKpVq2AymTLdNjo6Gu7u7rh06RJatWqFuXPnolmzZjleNzAwkNOBZEx9+8oMVFQUULq07mj+K/MfIABna0dRIHPnytby3r11R0IGt3nz5hw/v2TJEqxfvx5btmzJMjkBgPu/xYfd3NzQpUsXREREPDBBERnW8OHAkiVyLmrkSN3R5IrtTPFduCDl4/v2lSRFlE/h4eGYMWMG1q5di2LZ9A+Lj4/HzX8PN8bHx2PTpk3w9/e3ZphE5vXII1K3dO5c2QVtA2wnQYWGAsnJnNqjAhs2bBhu3ryJVq1aoX79+hg8eDAAICYmBu3atQMAXLx4EU2bNkW9evUQFBSE9u3bo02bNjrDJiq4ESOA06eBDRt0R5IrtrEGlZQEVKsmrwB++kl3NETZ4hoUGVpKClC9OuDnZ7RzUVnOs9vGCOqHH2SKj1vLiYjyz9lZdvRt2gQcPao7mgeyjQQ1Zw7g7Q2EhOiOhIjItr38MlC4MPDZZ7ojeSDjJ6h9+4D/+z9Zeypk/HCJiAytUiXg2WeBL78E4uN1R5Mj4z/jf/op8NBD3FpORGQuQ4cCN24Ay5frjiRHxk5Q164B33wj7YtLldIdDRGRfXj8cSAgAJg/H8h5o5xWxk5QS5YACQlSR4qIiMzDZAKGDAEOHAD27NEdTbaMm6DS0mQR7/HHgXr1dEdDRGRfevaUfnrz5+uOJFvGTVCbN0sPkyFDdEdCRGR/ihcH+vS52yHCgIyboD79FKhQQXabEBGR+Q0ZIhV6Fi7UHUmWjJmgzp4F1q0DBgwAXF11R0NEZJ9q1QKCg6VLeUqK7mgyMWaCCg2VnSWDBumOhIjIvg0dCpw7B6xfrzuSTIyXoJKSgC++ADp0kPp7RERkOR06AFWqGHKzhPES1KpVwKVL3BxBRGQN6fX5Nm8G/vlHdzT3MV6C+vRTwNOTdfeIiKxlwACpz/f557ojuY+xEtQffwA7dsjBXNbdIyKyjooVgWeeAcLCgDt3dEeTwVhZ4PPPZdde3766IyEiciyDBkl5uZUrdUeSwTgJKj4e+PproGtXoFw53dEQETmWp54CfHxky7lBGCdBrVgBxMVxazkRkQ4mkzz/7t4tyy0GYJwEFRoK+PoCTZvqjoSIyDH17i3LLAYZRRkjQR04AERESPY2ZdmanoiILK1cOeC554ClSw3RzNAYCWrBAqBIEeCll3RHQkTk2AYNkuWWb7/VHYkBEtStW8CyZUC3bkCZMrqjISJybI8/DtSpY4gzUfoT1PLlwM2b3BxBRGQEJhMweDCwdy+wf7/WUPQnqAULAH9/4NFHdUdCDmLChAkICAhA/fr1ERISgpiYmCxvFxYWBm9vb3h7eyMsLMzKURJp1KsXULSo9s0SJpVzP3rLNqvftw8IDATmzgWGDbPopYjSxcXFoWTJkgCAOXPm4MiRI/j8P9MZV69eRWBgICIjI2EymdCwYUPs27cPZR4wDZ3+NUQ2r29f4PvvgZgYoEQJS18ty91xekdQCxZIln7xRa1hkGNJT04AEB8fD1MWO0c3btyIVq1aoWzZsihTpgxatWqF8PBwa4ZJpNfgwbJH4JtvtIXgrO3KcXHywLt3B0qX1hYGOaZx48bhq6++QqlSpbB169ZMn4+OjkaVKlUy3vfw8EB0dHSW9xUaGorQ0FAAQGxsrGUCJrK2oCCgXj0ZSGjaI6BvBPXNN7LPnpsjyAKCg4Ph7++f6W3NmjUAgClTpuDcuXPo2bMn5s2bV6BrDRw4EJGRkYiMjESFChXMET6RfumVJfbvB37/XUsI+hLUF19Idg4K0hYC2a/Nmzfj8OHDmd46d+583+169uyJH374IdPXu7u749y5cxnvR0VFwd3d3eJxExnKCy/IMswXX2i5vJ4E9fvv8vbyy6wcQVZ37NixjP+vWbMGtWvXznSb1q1bY9OmTbh27RquXbuGTZs2oXXr1tYMk0i/0qWlskT6jJeV6UlQixZJ5YgePbRcnhzb22+/DX9/fwQEBGDTpk2YPXs2ACAyMhIDBgwAAJQtWxYTJkxAo0aN0KhRI7zzzjsoW7aszrCJ9BgwQPYMfP+91S9t/W3mt28DDz8MdOwo9Z6I7Ai3mZPdUQqoXRtwc5OGspZhkG3mP/wA3LghWZmIiIzNZJLn6507gb//tuqlrZ+gFi4EvLyAZs2sfmkiIsqHl14CnJ1lecaKrJugjh4Ftm8H+vfn5ggiIltRsSLQqRMQFgYkJVntstZNUIsWAU5O0hSLiIhsx4ABQGwssHat1S5pvQSVnAwsWQJ06ABUrmy1yxIRkRmEhABVqsgyjZVYL0GtXw9cusTNEUREtsjJCejXD9i0CThzxiqXtF6CWrhQtpe3aWO1SxIRkRn17Sv/fvmlVS5nnQR17hwQHi4PzllffVoiIiqAatVkqm/xYiA11eKXs06CWrIESEuT4SEREdmuAQNk0LFpk8UvZfkElZYmu/datgQ8PS1+OSIisqBOnYDy5a2yWcLy821pacCUKdy5R0RkD1xcgI8/BipVsvil9LZ8J7IzrMVHlC8GqcVHRESUC0xQRERkSExQRERkSExQRERkSExQRERkSExQRERkSExQRERkSExQRERkSExQRERkSExQRERkSA8qdUREeWAymcKVUmx6RmQGTFBERGRInOIjIiJDYoIiIiJDYoIiIiJDYoIiIiJDYoIiIiJD+n/sYrurR4h+RQAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<sympy.plotting.plot.Plot at 0x17f3be5acd0>"
     },
     "execution_count": 262,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#3.10\n",
    "\n",
    "#可以使用sympy.series函数\n",
    "import sympy as sym\n",
    "x = sym.symbols('x')\n",
    "formula = sym.cos(sym.sqrt(x**2+1))\n",
    "f1 = sym.series(formula,x,0,1).removeO()\n",
    "f3 = sym.series(formula,x,0,3).removeO()\n",
    "f5 = sym.series(formula,x,0,5).removeO()\n",
    "sym.pprint(f1)#1阶\n",
    "sym.pprint(f3)#3阶\n",
    "sym.pprint(f5)#5阶\n",
    "sym.plot((f1,(x,-3,3)),(f3,(x,-3,3)),(f5,(x,-3,3)),line_color = 'r')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.          0.          4.5         0.          8.91800933  0.85509858\n",
      "  12.99959308  2.75000597 16.35640054  5.7469784  18.62196582  9.63506456\n",
      "  19.74149169 13.99358057 20.03124502 18.48424233 19.97536056 22.9838953\n",
      "  20.08444862 27.48257286]\n",
      " [20.          0.         20.          3.         20.          6.\n",
      "  20.          9.         20.         12.         20.         15.\n",
      "  20.         18.         20.         21.         20.         24.\n",
      "  20.         27.        ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbiklEQVR4nO3de3RddZ338fc39zRJk7RJS29pKRRqGeiFTCkXEUWx8DCAN4QBQcexeGFGnoUPspzl5XE5z4ij6KOOlyosq1QBRYRRQBCdh9FCMb3RK7RAkrakado0TdM0t3O+zx9np4SQNGly9jlnJ5/XWlnnZO999v52n3M+/eW3Lz9zd0REJHqy0l2AiIiMjAJcRCSiFOAiIhGlABcRiSgFuIhIROWkcmMVFRU+Z86cVG5SRCTy1q1bd8DdK/tPT2mAz5kzh5qamlRuUkQk8sysbqDp6kIREYkoBbiISEQpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiIdp3uINvPPkirzS1JX3dCnARkRDVN7fznT/uYm/LsaSvWwEuIhKio509ABTlJ//CdwW4iEiI2oIALw4hwFN6LxQRkXHl8Ts5o/EIcIVa4CIikbJvMxNbtgNQnKcAFxGJlFg8MXB8UX520tetABcRCVHMnfycLHKykx+3CnARkRDF4h7KAUxQgIuIhCoW91AOYIICXEQkVHFXC1xEJJLUhSIiElGxeDhnoIACXEQkVOoDFxGJqJj6wEVEoimtLXAzm2VmfzKzbWa21cw+HUz/kpntNbONwc8VoVQoIhJRjhP38AJ8OGvtAW539/VmVgKsM7OngnnfdPevh1KZiEjE9V5GXxzSQcwhA9zdG4CG4PkRM9sOzAilGhGRMSTI78w4iGlmc4DFwNpg0q1m9oKZ3Wtm5YO8ZoWZ1ZhZTVNT0+iqFRGJkNdb4GkOcDMrBh4CbnP3VuD7wGnAIhIt9G8M9Dp3X+nu1e5eXVlZOfqKRUQi4vidCEO4lSwMM8DNLJdEeK92918DuHuju8fcPQ78CFgaSoUiIhF1vAVekL6zUAy4B9ju7nf3mT6tz2LvAbYkvzwRkeiKebhdKMNZ64XAh4DNZrYxmPY54HozWwQ4UAvcEkJ9IiKR9fpgDmkKcHf/M2ADzHos+eWIiIwdYY7GA7oSU0QkNPGQu1AU4CIiIeltgRfmqgUuIhIpsbiTbUbiXJDkU4CLiIQkFneys8IJb1CAi4iEJuZOlgJcRCR6ertQwqIAFxEJibpQREQiKq4AFxGJppg72eHltwJcRCQshuEhrl8BLiISErPXB3UIgwJcRCQkWWa4h5fgCnARkZCYQYj5rQAXEQlLlr1+Q6tQ1h/amkVExjkzHcQUEYkktcBFRCLKzNQHLiISRVmoBS4iEklqgYuIRJSpD1xEJJqy1AIXEYkmM4iHeCKhAlxEJCS9LfCwLqdXgIuIhKR3MJ6ekO5opQAXEQlJVpDg3bF4OOsPZa0iIkLvWA5dPQpwEZFI6W2Bpy3AzWyWmf3JzLaZ2VYz+3QwfZKZPWVmO4PH8lAqFBGJqN4+8M40tsB7gNvdfQGwDPiUmS0A7gSedvd5wNPB7yIiEkh7H7i7N7j7+uD5EWA7MAO4GlgVLLYKuCaUCkVEIqq3Bd6VCQcxzWwOsBhYC0x194Zg1j5g6iCvWWFmNWZW09TUNJpaRUQiJe194L3MrBh4CLjN3Vv7zvPEWeoDnujo7ivdvdrdqysrK0dVrIhIlPS2wNN6GqGZ5ZII79Xu/utgcqOZTQvmTwP2h1KhiEhEWZDgaTuIaYkK7gG2u/vdfWY9CtwcPL8ZeCT55YmIRFdWbx94SAGeM4xlLgQ+BGw2s43BtM8BXwUeNLOPAnXAtaFUKCISUUa4feBDBri7/5nXLyjq79LkliMiMnZkHe8D171QREQipbcPvCsWC2X9CnARkZCE3QeuABcRCcnrLXB1oYiIRIpa4CIiEWWZciWmiIicnJ54IriLC4ZzxvbJU4CLiISkszsR4LMnTQhl/QpwEZGQdHQnTh+cPVkBLiISKb33QJleVhjK+hXgIiIh6eiOkZ+TRW52OFGrABcRCUlHd5yCnOzQ1q8AFxEJSWdPjPzc8GJWAS4iEoLDx7rpiTsFuWqBi4hEyu7mdgAKctQCFxGJlLqDiQDPVwtcRCRa6pqPAqgLRUQkauoPtpOTZWTbYOPhjJ4CXEQkBPXN7aG2vkEBLiISirqD7aEewAQFuIhI0nX1xGk4fCzUA5igABcRSbo9h9qJe7gHMEEBLiKSdHXBOeD56kIREYmW4xfxqAUuIhItdQfbKcjNIjc7vFMIQQEuIpJ0dQfbqZo0AUMBLiISKfXNR6maVBT6dhTgIiJJ5O7UN7eHNoxaX0MGuJnda2b7zWxLn2lfMrO9ZrYx+Lki3DJFRKKh6UgnHd3xzAhw4CfA8gGmf9PdFwU/jyW3LBGRaOo9hbAqpJHo+xoywN39GaA59EpERMaA3tvIZkSAn8CtZvZC0MVSPthCZrbCzGrMrKapqWkUmxMRyXz1B4+SZTCzPHMD/PvAacAioAH4xmALuvtKd6929+rKysoRbk5EJBrqmtuZVlpIXshXYcIIA9zdG9095u5x4EfA0uSWJSISTak6AwVGGOBmNq3Pr+8Btgy2rIjIeFJ/MHUBnjPUAmb2C+ASoMLM9gBfBC4xs0WAA7XALeGVKCISDW2dPRw82sWsFBzAhGEEuLtfP8Dke0KoRUQk0uoOJsbBnJ2CqzBBV2KKiCRNfXAKYUb3gYuIyJs9vmUfBblZnFqhFriISGRs3nOYRze9xkcvOpWi/CF7p5NCAS4ikgR3PbGD8gm53PK201K2TQW4iMgoPfNSE3/edYBb3zGPiQW5KduuAlxEZBTiceerj+9gZnkhNy6rSum2FeAiIqPw6KbX2NbQymcuO5P8nHDHwOxPAS4iMkKdPTG+/uSLLJg2kasWTk/59hXgIiIjdN9z9ew5dIw7L59PVla4418ORAEuIjICrR3dfPePO7no9AouPiM9d1pVgIuIjMAP/9/LHGrv5rPL56etBgW4iMhJamzt4J4/v8pVC6dz9szStNWhABcROUnf+sNLxOLOZy47M611KMBFRE7Crv1tPPDX3dxw3myqUnTTqsEowEVETsLXntjBhLwc/ukdp6e7FAW4iMhwratr5sltjdxy8VwmF+enuxwFuIjIcLg7//bYDipL8vnoW09NdzmAAlxEZFj+sH0/NXWHuO2d85iQl5rbxQ5FAS4iMoSeWJy7ntjB3IoiPlg9K93lHKcAFxEZwkPr97Brfxt3LD+TnOzMic3MqUREJAMd64px91MvsaSqjHefdUq6y3kDBbiIyAnc+5dXaWzt5M7L34JZ6m9YdSIKcBGRQRw62sUP/utl3vmWKSw9dVK6y3kTBbiIyCC++6ddHO3q4Y403rDqRBTgIiID2N3czs+ereP9587kjKkl6S5nQApwEZEB3P3US5jB/3zXGekuZVAKcBGRfra+dpjfbNzLRy48lWmlhekuZ1BDBriZ3Wtm+81sS59pk8zsKTPbGTyWh1umiEjq3PXEi5QW5vKJS05LdyknNJwW+E+A5f2m3Qk87e7zgKeD30VEIu8vuw7wzEtN3Pr20yktzE13OSc0ZIC7+zNAc7/JVwOrguergGuSW5aISOrF485XH9/BjLJCblw2O93lDGmkfeBT3b0heL4PmDrYgma2wsxqzKymqalphJsTEQnf7zY3sHnvYW6/7AwKcrPTXc6QRn0Q090d8BPMX+nu1e5eXVmZnpGbRUSG0tUT599//yLzTynh6kUz0l3OsIw0wBvNbBpA8Lg/eSWJiKRWR3eMz/xyE/XN7Xz28vlkZ2XWJfODGWmAPwrcHDy/GXgkOeWIiKTWay3H+MAPnuU/X3iNO5afySVnRKenYMi7kpvZL4BLgAoz2wN8Efgq8KCZfRSoA64Ns0gRkTDU1Dbz8fvW09Ed48c3VXPpWwY9nJeRhgxwd79+kFmXJrkWEZGUuf/5ej7/yBZmlBVy/4rzOH1KZl4ufyKZMS6QiEiKdMfifOW321j1bB1vnVfBd69fQumEzD7fezAKcBEZN5qPdvGp1et59pWDrLh4Lne8O7NG2DlZCnARGRd27GvlH1fVsP9IJ3dfu5D3LpmZ7pJGTQEuImPe45sbuP2XmygpyOHBW85n0ayydJeUFApwERmz4nHnW0/v5NtP72RxVRk/vPFcpkwsSHdZSaMAF5Exqa2zh9sf3Mjvtzby/nNn8pVr/iYSl8efDAW4iIw59Qfb+dhPa9jV1MYXrlzARy6ck3EDEieDAlxExpS/7DrAp36+HndY9ZGlXDSvIt0lhUYBLiJjgrvzkzW1fOV32zmtsogf3VTN7MlF6S4rVApwEYm8zp4Yn//NFh6s2cO7Fkzlmx9cRHH+2I+3sf8vFJExbX9rB7fct44N9S3886XzuO3SeWRF5G6Co6UAF5HI2rS7hVt+to7Dx7r53g1LuOLsaekuKaUU4CISSQ9v2MNnH9rMlJJ8HvrEBSyYPjHdJaWcAlxEIiUWd+56Ygcrn3mFZXMn8b0bzmVSUV66y0oLBbiIRMbh9m7+6f4NPPNSEzedP5vPX7mA3AjfjGq0FOAiEgm79h/hYz9dx55D7fzbe8/m+qVV6S4p7RTgIpLxnt7eyKfv30hBbhY//9gy/nbOpHSXlBEU4CKSsdyd7/3Xy3z9yRc5a/pEfvihamaUFaa7rIyhABeRjHSsK8b/+tUmfvtCA3+3cDpfe985FOaNrZtRjZYCXEQyzt6WY6z4aQ3bGlr57PL5fPxtc8fkzahGSwEuIhnl+Veb+cR96+jqiXPPzdW8Y360RopPJQW4iGSM1Wvr+OIjW6maNIGVN1Vz+pTidJeU0RTgIpJ23bE4//s/t3Lfc/W87YxKvn39YkoLozlSfCopwEUkrQ60dfLJ1et5/tVmbrl4Lncsn0/2OLkZ1WgpwEUkLV5qPMKqNbU8vGEvsbjzrQ8u4prFM9JdVqQowEUkZXpicf6wfT+r1tTy7CsHycvJ4uqF0/nYxXM5Y2pJusuLHAW4iISu+WgX9/+1ntXP1bO35Rgzygr57PL5fPBvZ43bG1Elw6gC3MxqgSNADOhx9+pkFCUiY8PmPYdZ9Wwtj256ja6eOBecNpkv/N0CLp0/hZxxfBOqZElGC/zt7n4gCesRkTGgqyfO41saWLWmlvX1LUzIy+ba6pncfP4c5qmbJKnUhSIiSdHY2sHqtfX8fG09B9o6ObWiiC9cuYD3V89kYoFOCQzDaAPcgSfNzIEfuvvK/guY2QpgBUBVlW7/KDKWuDs1dYdYtaaWJ7bsI+bO28+cwk3nz+bieZXjZmzKdBltgF/k7nvNbArwlJntcPdn+i4QhPpKgOrqah/l9kQkA3R0x3hk415WraljW0MrJQU5fPiCOdy4bDZzKorSXd64MaoAd/e9weN+M3sYWAo8c+JXiUhU7W5u577n6nigZjct7d2cObWE//Oes7lm8XQm5KlHNtVGvMfNrAjIcvcjwfPLgC8nrTIRyQjuzl92HeQna2p5ekcjWWZctmAqN18wh/NOnaS7BKbRaP7LnAo8HLx5OcDP3f2JpFQlImnX1tnDr9fvYdWaWl5uOsqkojw+eclp3HDebKZrUIWMMOIAd/dXgIVJrEVEMsDLTW387Nk6frVuD22dPZwzs5RvfGAh/+OcaRTkakCFTKJOKxEhFnf+tGM/q56t5b93HiA327jynOncdP5sFleVp7s8GYQCXGQca2nv4sGa3fzsuTp2Nx9j6sR8bn/XGVy3tIrKkvx0lydDUICLjEPbG1pZtaaW32zcS0d3nKVzJnHn8rdw2VlTydUl7pGhABcZJ7pjcZ7c2siqNbU8X9tMQW4W1yyawU3nz2HB9InpLk9GQAEuMsY1HenkF8/Xs3ptHY2tncwsL+RzV8zn2upZlE3QnQCjTAEuMkZtqE9c4v67zQ10x5y3zqvgX685m7fPn6IRb8YIBbjIGNLRHeN3LzTw02dr2bTnMMX5Odxw3mxuXDZbAwSPQQpwkTHgtZZjrF5bx/3P7+bg0S5Oqyziy1efxXsWz6BEdwIcsxTgIhHl7qx9tZlVa2p5clsjcXcunT+VD18whwtPn6xL3McBBbhIhLR2dLNpdwsb6lt4bHMDO/YdobQwl3+86FRuXDabWZMmpLtESSEFuEiGisedXU1tbKg/xPq6FjbsPsTO/W24gxmcPaOUu953NlctnEFhni5xH48U4CIZoqW9iw1B63pD/SE21rdwpLMHgNLCXBZXlXHlOdNZXFXGwlllGuVGFOAi6dATi/NSYxsbdr/eun6l6SgAWQZnnjKRqxZNZ3FVOUuqyji1okh92vImCnCRFDjY1smG+hbW1x9iQ30Lm/a00N4VA2ByUR6Lq8p535KZidb1zDKK8vXVlKHpUyKSZN2xODsajgRhfYgNu1uoO9gOQE6WsWD6RD5w7sygdV3OrEmFal3LiCjARUZpf2vH8Zb1hvoWXtjbQkd3HIApJfksqSrn75dWsWR2OX8zvVQHHCVpFOAiJ6GzJ8a211pZHxxo3FDfwt6WYwDkZWdx1oyJ/P3S2SyuKmPJ7HKmlxaodS2hUYCLDMLdaTj8eut6ff0htu5tpSuWaF1PLy1g8exyPnLhHJbMLues6RPJz1HrWlJHAS4S6OiOsXnv4Tecd93Y2glAfk4W58ws5cMXzmFJVRmLZpVzSmlBmiuW8U4BLuOSu7O7+Rgbdr/eut72Wis9cQegatIEls2dzOJZia6Q+adMJC9HAx1IZlGAy7hwtLOHF/YcPn7e9cbdhzjQ1gVAYW42C2eV8rGL57KkqpxFs8o0nJhEggJcxhx359UDR99w3vWOfa0EjWvmVhRx8RmVLKkqZ3FVGWdOLSFHw4hJBCnAJfKOdHSzaffhN5x33dLeDUBxfg6LZpVx69tPZ3HQui4v0ig0MjYowCUjxeJO89EuDrR10nSkc4DHruO/Hzzadfx186YUc9mCqUHrupzTpxRr9BkZsxTgkjLxuNNyrPtNYdzU+9gnmJuPdh7v8uirIDeLypJ8KorzmT15AufOKWfaxAIWzkrc4Km0UDd4kvFDAS6j4u60Huuhqa2DpiNdNLV1ciAI5eOPQUAfbOs6fpZHX3k5WVQW51NRks+MsgIWzSqlojj/eFD3fSzKy9aFMSIBBbi8ibvT1tnzpq6K/q3mA8H83gtb+srJsuOhO6WkgAXTJg4YyBXF+UwsyFEoi4zAqALczJYD/xfIBn7s7l9NSlUSiqOdPQME8cAB3dnz5lDOzjImF+UdD995U0qCEM6jsiSfyj6hXFqYS5b6nkVCNeIAN7Ns4D+AdwF7gL+a2aPuvi1ZxcnQOrpjb2gRJx67gi6NN7age29f2pcZbwjluRVFVARhXFGSR2VxQfCYT/mEPIWySAYZTQt8KbDL3V8BMLP7gauBpAf4d57eyaObXkv2aiOtOxbnYFvX8RFb+iufkHu8Nby4qmyArotEq3nShDydAy0SUaMJ8BnA7j6/7wHO67+Qma0AVgBUVVWNaEOVJfnMm1o8oteOVdlZWUwuemPXRW9ATy7OI1ehLJJ+p5wd6upDP4jp7iuBlQDV1dUDnBg2tOuWVnHd0pGFv4hI2lwe7mHB0TTT9gKz+vw+M5gmIiIpMJoA/yswz8xONbM84Drg0eSUJSIiQxlxF4q795jZrcDvSZxGeK+7b01aZSIickKj6gN398eAx5JUi4iInASdqiAiElEKcBGRiFKAi4hElAJcRCSizH1E19aMbGNmTUDdCF9eARxIYjnJorpOjuo6OZlYVybWBGO7rtnuXtl/YkoDfDTMrMbdq9NdR3+q6+SorpOTiXVlYk0wPutSF4qISEQpwEVEIipKAb4y3QUMQnWdHNV1cjKxrkysCcZhXZHpAxcRkTeKUgtcRET6UICLiERUxgW4mS03sxfNbJeZ3TnA/HwzeyCYv9bM5qSgpllm9icz22ZmW83s0wMsc4mZHTazjcHPF8KuK9hurZltDrZZM8B8M7NvB/vrBTNbkoKazuyzHzaaWauZ3dZvmZTsLzO718z2m9mWPtMmmdlTZrYzeCwf5LU3B8vsNLObU1DXv5vZjuB9etjMygZ57Qnf8yTX9CUz29vnfbpikNee8HsbQl0P9Kmp1sw2DvLaUPZVsO4BcyGlny93z5gfErelfRmYC+QBm4AF/Zb5JPCD4Pl1wAMpqGsasCR4XgK8NEBdlwC/TcM+qwUqTjD/CuBxwIBlwNo0vKf7SFyIkPL9BVwMLAG29Jn2NeDO4PmdwF0DvG4S8ErwWB48Lw+5rsuAnOD5XQPVNZz3PMk1fQn4zDDe4xN+b5NdV7/53wC+kMp9Fax7wFxI5ecr01rgxwdKdvcuoHeg5L6uBlYFz38FXGpmoQ6V7u4N7r4+eH4E2E5iTNAouBr4qSc8B5SZ2bQUbv9S4GV3H+kVuKPi7s8Azf0m9/0MrQKuGeCl7waecvdmdz8EPAUsD7Mud3/S3XtHqX6OxChXKTPIvhqO4XxvQ6kr+O5fC/wiWdsbrhPkQso+X5kW4AMNlNw/KI8vE3zYDwOTU1IdEHTZLAbWDjD7fDPbZGaPm9lZKSrJgSfNbJ0lBpDubzj7NEzXMfiXKx37C2CquzcEz/cBUwdYJt377R9I/OU0kKHe82S7NejWuXeQ7oB07qu3Ao3uvnOQ+SnZV/1yIWWfr0wL8IxmZsXAQ8Bt7t7ab/Z6Et0EC4HvAL9JUVkXufsS4HLgU2Z2cYq2OyRLDLV3FfDLAWana3+9gSf+ns2oc2nN7F+AHmD1IIuk8j3/PnAasAhoINFdkUmu58St79D31YlyIezPV6YF+HAGSj6+jJnlAKXAwbALM7NcEm/Sanf/df/57t7q7m3B88eAXDOrCLsud98bPO4HHibx52xf6Rx8+nJgvbs39p+Rrv0VaOztRgoe9w+wTFr2m5l9GLgSuCH48r/JMN7zpHH3RnePuXsc+NEg20rXvsoB3gs8MNgyYe+rQXIhZZ+vTAvw4QyU/CjQe8T2/cAfB/ugJ0vQz3YPsN3d7x5kmVN6++LNbCmJfRvqfyxmVmRmJb3PSRwE29JvsUeBmyxhGXC4z593YRu0dZSO/dVH38/QzcAjAyzze+AyMysPug0uC6aFxsyWA3cAV7l7+yDLDOc9T2ZNfY+XvGeQbaVrgPN3Ajvcfc9AM8PeVyfIhdR9vsI4OjvKI7tXkDia+zLwL8G0L5P4UAMUkPiTfBfwPDA3BTVdROLPoBeAjcHPFcDHgY8Hy9wKbCVxBP454IIU1DU32N6mYNu9+6tvXQb8R7A/NwPVKXofi0gEcmmfaSnfXyT+A2kAukn0M36UxDGTp4GdwB+AScGy1cCP+7z2H4LP2S7gIymoaxeJftHez1jv2VbTgcdO9J6HWNPPgs/NCySCaVr/moLf3/S9DbOuYPpPej9PfZZNyb4K1j9YLqTs86VL6UVEIirTulBERGSYFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYj6//Q+8cmsIkaFAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#3.11\n",
    "\n",
    "#最终计算结果为9秒钟左右\n",
    "import numpy as num\n",
    "import numpy.linalg as ng\n",
    "import matplotlib.pyplot as plt\n",
    "N = 2 #需要更新位置的即是兔子和猎狗两个对象\n",
    "vr = 3#兔子的速度为3m/s\n",
    "vd = 4.5 #猎狗的速度为4.5m/s\n",
    "times = 9#设定初步的追逐时间\n",
    "divs = 10\n",
    "T = np.linspace(0,times,divs)\n",
    "# print(T)#可以简要查看一下位置的更新迭代时间\n",
    "dt = T[1] -T[0]\n",
    "xy = np.array([[0,0],[20,0]])#初始位置，猎狗与兔子\n",
    "Txy = xy#Txy用来存储所有的xy坐标\n",
    "xyn = np.empty((2,2))#空矩阵用来接收更新的兔子和猎狗的位置\n",
    "for n in range (1,len(T)):#每个时间点一次循环\n",
    "    for i in [0,1]:#只有两个待更行对象\n",
    "        j = (i+1)%2\n",
    "        dxy = xy[j]-xy[i]#两者的相对位置\n",
    "        dd = dxy/ng.norm(dxy)#单位化向量，相当于给与了运动的方向\n",
    "        if i ==0:\n",
    "            xyn[i] = xy[i] + vd*dt*dd #猎狗的位置更新\n",
    "        else:\n",
    "            xyn[i] = [20,xy[i][1]+vr*dt]#兔子的位置更新\n",
    "    Txy = np.c_[Txy,xyn]#矩阵存储当前的坐标\n",
    "    xy = xyn.copy()#更新时刻点的xy坐标\n",
    "print(Txy)\n",
    "for i in range(N):\n",
    "    plt.plot(Txy[i,::2],Txy[i,1::2])\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "一重积分计算结果为： (0.8328534212790222, 1.2431818039574587e-08)\n",
      "二重积分的值为： 21.40714285713916\n",
      "三重积分计算结果为： (67.02064327658225, 7.440786129085081e-13)\n"
     ]
    }
   ],
   "source": [
    "#3.12\n",
    "\n",
    "#（1）\n",
    "import numpy as np\n",
    "from scipy.integrate import quad,dblquad\n",
    "f = lambda x: np.exp(-x)*np.sin( np.sqrt(x**2+2))\n",
    "print('一重积分计算结果为：', quad(f,0,np.inf))#返回值有两个，第一个值是积分的值，第二个值是对积分值的绝对误差估计。\n",
    "#（2）对于2问中，可以将积分划分为（0,1）及（1,4）区域\n",
    "f1 = lambda y,x:x**2+2*y**2\n",
    "bd1 = lambda x:np.sqrt(x)\n",
    "bd2 = lambda x:x-2\n",
    "value1,err1 = dblquad(f1,0,1,lambda x: -bd1(x),bd1)#输出值和误差\n",
    "value2,err2= dblquad(f1,1,4,lambda x:bd2(x),bd1)\n",
    "print(\"二重积分的值为：\",value1+value2)\n",
    "#（3）三重积分\n",
    "# x**2+y**2= p**2 = z，说明曲面半径为np.sqrt（z）,则曲面面积为np.pi*z\n",
    "print('三重积分计算结果为：', quad(lambda z:z*np.pi*z,0,4))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2*exp(-x)*sin(x) + 2*exp(-x)*cos(x)\n",
      "⎧        5⋅π │      ⎫   ⎧        π │      ⎫\n",
      "⎨2⋅n⋅π + ─── │ n ∊ ℤ⎬ ∪ ⎨2⋅n⋅π + ─ │ n ∊ ℤ⎬\n",
      "⎩         4  │      ⎭   ⎩        4 │      ⎭\n",
      "极大值为：0.644793883889669\n",
      "极小值为：-0.027864070195388408\n"
     ]
    }
   ],
   "source": [
    "#3.13\n",
    "\n",
    "#求极值点\n",
    "from scipy.optimize import fminbound\n",
    "x = sym.symbols('x')\n",
    "expr1 = 2*sym.exp(-x)*sym.sin(x)\n",
    "expr1_x = sym.diff(expr1,x)\n",
    "print(expr1_x)#输出一阶导\n",
    "sym.pprint(sym.solveset(expr1_x,x))#输出一阶导为0的点，可见三角函数的解有无穷多个\n",
    "#当n为0时，一阶导为0的点为sym.pi/4,5*sym.pi/4\n",
    "f= lambda x:2*np.exp(-x)*np.sin(x)\n",
    "print('极大值为：{}'.format(f(np.pi/4)))\n",
    "print('极小值为：{}'.format(f(5*np.pi/4)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.27433388230814\n",
      "28.27433388230814\n",
      "56.54866776461628\n",
      "554176.9440932396\n"
     ]
    }
   ],
   "source": [
    "#3.14\n",
    "\n",
    "#旋转一周形成葫芦一样的柱体\n",
    "#（1）求体积，将图像分割为-2<=y<=1,1<y<=4两部分\n",
    "from scipy.integrate import quad\n",
    "func1 = lambda y:np.pi*(4-y**2)\n",
    "func2 = lambda y:np.pi*(4-(y-2)**2)\n",
    "val1,er1 = quad(func1,-2,1)#下半部分旋转的体积\n",
    "print(val1)\n",
    "val2,er2 = quad(func2,1,4)#上半部分旋转的体积\n",
    "print(val2)\n",
    "print(val1+val2)#最终的合并体积  #可以整体使用sympy求符号解\n",
    "#（1）求做功(功 = 力*距离，水的质量 = 水的体积*水的密度，力F = mg，可以先求单位体积的做功，再求整体的微分)\n",
    "p = 1*10**3 #水的密度\n",
    "g = 9.8 #重力加速度\n",
    "func3 = lambda y: p*g*(2-y)*np.pi*(4-y**2)\n",
    "func4 = lambda y: p*g*(2-y)*np.pi*(4-(y-2)**2)\n",
    "val3,err3 = quad(func3,-2,1) #抽取下班部分水所做的功\n",
    "val4,err4 = quad(func4,1,4) #抽取上半部分水做的功\n",
    "print(val3+val4)#最终计算出来的做功"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "阻力系数为： 7000.0\n",
      "安全距离为: 685.7142857142857\n"
     ]
    }
   ],
   "source": [
    "#3.15\n",
    "\n",
    "#参考文档：https://www.docin.com/p-1853085899.html\n",
    "#可以不进行单位的换算，直接计算\n",
    "#（1）\n",
    "'''\n",
    "飞机的速度为V，则减速伞的阻力f = -kv\n",
    "牛顿第二定律F（合） = ma,应用牛顿第二定律可得到：m*dv/dt = -kv\n",
    "为得到速度与距离之间的关系将dv/dt改写为dv/dx*dx/dt = v*dv/dx\n",
    "将此式子代入上牛顿第二定律推导式，则dv = -k/m*dx，\n",
    "对左右两边同时积分(左边vo--vt，右边想x = x（t）),求出K = m*（v0-vt）/x(t)\n",
    "'''\n",
    "m = 5000 #飞机重量\n",
    "v0 = 800 #初速度为800\n",
    "x = 500 #飞机滑行的距离\n",
    "vt =100 #飞机最后的末速度\n",
    "k = m *(v0-vt)/ x\n",
    "print('阻力系数为：',k)\n",
    "#(2)\n",
    "'''\n",
    "即相同的阻力系数\n",
    "由于m*dv/dt = -kv，则dv/v =- k/m*dt,同样应用上问中的定义域积分方法\n",
    "求解出v(t),再利用dx/dt = v(t)，再次积分，求出x（t）\n",
    "'''\n",
    "m2 = 8000\n",
    "v2 = 600\n",
    "f = lambda t :m2*v2/k*(1-np.exp(-k/m*t))#该方程小于等于m2*v2/k\n",
    "print('安全距离为:',m2*v2/k)#685<1200,所以能够安全着陆"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "局部极小点为：[1.31759894 1.73607197]，极小值为：-0.00016724663812221996\n"
     ]
    }
   ],
   "source": [
    "#3.16\n",
    "\n",
    "#求二元函数的局部极小点\n",
    "from scipy.optimize import minimize\n",
    "f = lambda x: 100*(x[1]-x[0]**2)**2+(1-np.sin(x[0]))**2*np.cos(x[1])\n",
    "x0 = minimize(f,[2.0,2.0],)\n",
    "print('局部极小点为：{}，极小值为：{}'.format(x0.x,x0.fun))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}